program test_tbrconv

    use iasing_osse
    use random_numbers_mix
    use netcdf
    
    implicit none

    integer :: ionc, ncid
    character(len=128) ::  arg(5)
    integer :: iatm, natm, ilast
    integer :: open_mode

    type(r4a_ids) :: db_ids
    real(JPR4), dimension(NIASI) :: wia, ria, tia
    real(JPR4), dimension(NIANG) :: wng, rng, tng
    real(JPR4), dimension(NIANG) :: rnga, rngb, tnga, tngb
    logical :: ngb_existence

    type(randomNumberStream) :: randstream
    real(JPR8) :: zrand(NIANG)
    integer(JPI4) :: kseed1, kseed2, kseedT
    integer(8) :: t8, rate8, cmax8


    include "r4a_initializer.interface"
    include "rad2brt_f.interface"
    include "brt2rad_f.interface"
    include "checknc.interface"

    ! Get Last simulation
    open(1,file='last',action='read',status='old')
    read(1,*) ilast
    close(1)
    if (ilast.eq.0) then
        print*, "bad last value. Stop"
        stop -1
    else
        print*, "ilast: ", ilast 
        ilast = ilast + 1
        print*, "ilast: ", ilast 
    endif
    !stop


    ! Initialisation of variables and modes
    open_mode = NF90_WRITE
    kseed1 =     524287
    kseed2 = 2147483647

    call get_command_argument(1,arg(1))
    print*, trim(arg(1))

    ionc = nf90_open(arg(1), open_mode, ncid)
    if (ionc /= NF90_NOERR) then
        print*, "Netcdf error: ", nf90_strerror(ionc)
        print*
        stop
    endif

    ! Initialise databank
    db_ids = r4a_initializer(ncid)

    ! Check if ngb field exists or not in the databank
    ngb_existence = .false.
    if (db_ids%rbt%gid .ge. 0) then
        ! Case where 'rbt' group exists
        if (db_ids%rbt%ngb .gt. 0) ngb_existence = .true.
    else
        print*, "Databank error: 'Radiances' group does not exists..."
        print*, "Aborting execution"
        stop -1
    endif
    ! Create ngb variable
    if (open_mode .eq. NF90_WRITE .and. .not. ngb_existence) then
        call checknc(nf90_redef(ncid))
        call checknc(nf90_def_var(db_ids%rbt%gid, "ngb", NF90_FLOAT, &
                          (/db_ids%wng_did, db_ids%atm_did/), db_ids%rbt%ngb))
        call checknc(nf90_enddef(ncid))
    endif
    



    ! Inquire length of sims
    ionc = nf90_inquire_dimension(ncid, db_ids%atm_did, len=natm)

    ! Get wavenumber arrays
    ionc = nf90_get_var(db_ids%dim%gid, db_ids%dim%wia, wia)
    ionc = nf90_get_var(db_ids%dim%gid, db_ids%dim%wng, wng)


    do iatm = ilast, natm

        ! Get radiances
        ionc = nf90_get_var(db_ids%rbt%gid, db_ids%rbt%ria, ria, &
            start=(/1,iatm/), count=(/NIASI,1/) )
        ! Cycle if simulation not run
        if (ria(1) <= 0) cycle

        ionc = nf90_get_var(db_ids%rbt%gid, db_ids%rbt%rng, rng, &
            start=(/1,iatm/), count=(/NIANG,1/) )
        
        ! Convert radiances to brightness temperature
        if (ria(1).gt.1) then
            ! Change units from mWatts to Watts
            ria = ria/1E3
            rng = rng/1E3
        endif
        tia = rad2brt_f(NIASI, wia, ria)
        tng = rad2brt_f(NIANG, wng, rng)

        !print*, " -- ria: ", ria(1)/1E3
        !print*, " -- tia: ", tia(1)

        ! Initialise random seeds
        kseedt = modulo(int(kseed1,JPI8) + int(kseed2,JPI8), &
            int(huge(kseedt), JPI8))
        kseed2 = kseed1
        kseed1 = kseedt
        call initialize_random_numbers(kseed1, randstream)

        ! Get random seeds and modify brightness temperatures
        call gaussian_distribution(zrand(1:NIASI), randstream)
        tia  = tia + real(zrand(1:NIASI),4)*NO_IASI
        call gaussian_distribution(zrand(1:NIANG), randstream)
        tnga = tng + real(zrand,4)*NO_IANGA
        call gaussian_distribution(zrand(1:NIANG), randstream)
        tngb = tng + real(zrand,4)*NO_IANGB
        
        ! Return to radiance space 
        !print*, " -- tia: ", tia(1)
        ria  = brt2rad_f(NIASI, wia, tia)
        !print*, " -- ria: ", ria(1)
        rnga = brt2rad_f(NIANG, wng, tnga)
        rngb = brt2rad_f(NIANG, wng, tngb)

        if (open_mode .eq. NF90_WRITE) then
            ionc = nf90_put_var(db_ids%rbt%gid, db_ids%rbt%ria, ria, &
                start=(/1,iatm/), count=(/NIASI,1/))
            ionc = nf90_put_var(db_ids%rbt%gid, db_ids%rbt%rng, rnga, &
                start=(/1,iatm/), count=(/NIANG,1/))
            ionc = nf90_put_var(db_ids%rbt%gid, db_ids%rbt%ngb, rngb, &
                start=(/1,iatm/), count=(/NIANG,1/))
        endif
        print*, iatm

        !if (iatm==100000) exit
    enddo



end program test_tbrconv
